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ABSTRACT 


This thesis is concerned with the numerical 
approximation of initial value problems in differen- 
tial equations. 

A detailed development of Nordsieck's method 
is presented. The errors introduced by stepsize 
changes are then e€xamined and a procedure to reduce 
these errors is introduced. Theoretical and practical 
considerations showed that the procedure is effective . 


in reducing errors. 
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CHAPTER I 
INTRODUCTION 
dNewob) ecculveroOL this thesiss#is to study 


Nordsieck's method for solving the differential 


equation 


yee eesy) & (1.1.1) 


given y(t) me A 
Nordsieck's method consists of a predictor 


equation of the form 
ce ORRE oocat AzGeoy 


and a corrector equation of the form 


(mils) eee) (m) 2. 
eS ahr, Pini hare Sey ae aI ae aN i ar 
Cie tees) 
where 
at 2 
x, = col(y, shy} ae Bena tidy auc 
(Ge Pe argc) 


a) 


Here, vi is the j-th derivative of a polynomial that 
approximates the solution about a point t.- The matrix 


Q and vector goare chosen so that the eruncatwon error 


(m) 


is low and the method is stable. The quantities Ey 


(2.1.2) 
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are scalars that measure how well the differential equa- 
tion is satisfied. 

One of the principal advantages of Nordsieck's 
method is that estimates of local errors are easily 
obtainable. In fact, when a q-th order method has been 
usedstoumakesassteptthenmthe erroraforithe (q-1)i<th 
order method is immediately available since all deriva- 
tives from the first to the q-th are available; and by 
taking differences estimates of the errors for the g-th 
and (q+tl)-th order methods can easily be computed. 

Because local estimates are readily available it 
is possible to choose both order and stepsize to mini- 
mize, at least locally, the cost of solving differential 
equations for a given tolerance. 

An implementation of Nordsieck's method with 
facilities for changing the stepsize and order has been 
written and is included in the Appendix. While experi- 
menting with this program we observed that large errors 
are introduced into the numerical solution when using 
methods of order greater than 5. Similar observations 
are also reported by Krough [15]. Hence, we developed 
and included in the program a procedure to substan- 
tially lower the stepsize changing errors when the step- 


size decreases. 
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TpeChapterelithie predictor ‘and ‘corrector 
equations of Nordsieck's method are developed using 
matrix notation. Matrix notation simplifies the 
analysis considerably, it illustrates the method's 
relation to multistep methods more readily, and it 
leads quite easily to a more general class of inte- 
gration methods called multivalue methods by Gear [7]. 

In Chapter III stepsize changing is discussed. 
imeparticular, Gitrerent methods of Obtaining local 
truncation error estimates are considered. Also, we 
present the development of a technique to correct the 
solution when decreasing the stepsize. 

In Chapter IV the implementation of Nordsieck's 
method is discussed. Chapter V contains the test pro- 


blems. 
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NORDSIECK'S METHOD 


ea MNO CACLON 


het jy be Che true) solution and Y be the appro- 


ximatesolutron=of (1.1.1). We use the abbreviations 
Need ARE Oh ee 
oe = Y(t.) I 

and yo = E(t iY.) ' 
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WHheleerc sot 
n 0 


Dal 
the stepsize. 
Vectors are denoted by underlined lower case 


letters. The symbol ||.|| is used for the Euclidean 


VecLOGMnOrl. Ome xanpLe, 


x") (t) = A x(t) 


is a constant coefficient linear system, x(t) 5 is the 
i-th component of the solution vector xO) 7 and 
I fx(t) ||? = F(x (t),)° is the square of the Euclidean 


norm. Moreover, if Xj1 Xor X30 ceer x is a sequence 
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of vectors,. then x4 denotes the i-th component of the 
v 
n-th vector in the sequence. 
For convenience zero indexing is wsed for vectors 


and matrices. 
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Tosdeveloperie=predictor equation (1.1.2) OL 
Nordsieck's method we consider the truncated Taylor 


series 
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assuming that the approximate solution Y has derivatives 


Of SUEFICiently high orders over the interval [t_,t 
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tion by using the vectors 
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Relationships between x, and x 
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Xx, and Zo, can be obtained by expanding the components 


Xe] and 254] 17 truncated Taylor series about mae 
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the use of "(2.2.6)2 Furthermore, if Zz. Constscts "OL the 
derivatives of a q-th degree polynomial evaluated at 
t = toe then z consists of the derivatives of the 


n+1 


same polynomial evaluated at t =t Recalling that 


n+1° 


the derivatives.of a polynomial evaluated at a.point 
uniquely define the polynomial, we see that Zz, and 241 
are simply different representations of the same poly- 


nomial defined by 


a ” 
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4=0 aes a 

Since the sequence Zor Zyr cee Simply gives 


different representations of the polynomial defined 


by Zz Le should be clear=that*the use of the predictor 


0. 
alone cannot yield good approximations to the solution 
Of a differential equation. “Example 2.iTirrrlustrates 
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Using this polynomial one can obtain the second column 
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An alternate method of coming to the same con- 
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2.4 Parameter Definitions 


In section 2.3, immediately following equation 
(2.3.16), we choose & = toler: Unfortunately, 
using Chis) choice for 2 in the corrector equation 
(2.3.1) produces an unstable method. However, there 
SxS eevectOLos swiichemakes the method (2.3.1 )as cable. 
One way to determine such an % is to show that the 
corrector equation (2.3.1) when iterated to convergence 
is equivalent to a multistep method and then to impose 
the stability condition of multistep methods on (2.3.1). 
We use Osborne's [17] proof to show that the stability 
O£ the multivalue method (2.3.5) depends only on 
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To determine a particular k we use the following 
two lemmas, whose proofs are deferred to the end of 
thasesection, to show that (2.4.1) is equivalent to a 


multistep method whose parameters depend on k. 
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linear combination of all the Y values present in the 
equation. Gn and I, are linear combinations of f 
values. Moreover, if id can be made identically zero, 


then the resulting equation 
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is a multistep method. One way to make Jo, equal zero 
is to choose the WR to be the coefficients of the 
characteristic polynomial of A. With this choice the 


Hamilton-Caley theorem asserts 
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Thus, (2.4.13) becomes identically zero and (2.4.9) 
becomes the multistep method (2.4.14). 
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a (E)Y¥, = hp (E)f, (26142, 18) 


where q and g are polynomials. Moreover, we recall that 
given both the q-polynomial of a multistep method and 
its order, the g polynomial is determined. Furthermore, 
aed CLStep method) is said, to; be stable fi. allaroots.of 
the a-polynomial are less than or equal to one in magni- 
tude and if those roots that are equal to one are simple. 
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polynomial w(z) defined.by, (2.4:17).; Now, recali that 
Wize ws thescharacteristichequati Ons ofs An—Bi- k dar 
Moreover, recall that A has been obtained from Onn Sloe Q) 
by deleting the first two rows and columns. With these 
remarks we use the following theorem stated and proved 


by.Gear, [5]..to assert. the, existence, of, a, unique k. 


Theorem 
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numbers {r,} a column vector 2 exists such that the 


matrix 


M 
= iG 
rr See Sato 


has i eigenvalues equal to one and q-it+l eigenvalues 
equal to the set {A;}- The eigenvalues are independent 
of the first 1 elements of 2 and’ tnigquely detemmine gthe 
TAS tel pies Fements sot, 2 ¢ 

In the remainder of this study we assume that 
eS OO, eb=Hl,2;....fq-L, where tA; J is the set of roots 


of y(z). Consequently, 


atthy Se) Sai) ere (2.4.20) 


aAT* = 0, n. = 1, and Hes Op for inl, 2). Sfqal se More= 


g 
OVeEry wile thvs echorce of {A; } Nordsieck's method is 


stable. Furthermore, in the next section it is shown 
that Qo can be chosen so that the local truncation error 


soe Hence, from 


of the (q+l) value method is O(h 
previously made remarks it follows that Nordsieck's 
method is equivalent to the Adams-Moulton method. 


The proofs of the lemmas complete this section. 


Proof of Lemma 2.4.1 


Write equations (2.4.la) and (2.4.1b) in the form 
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el ane 


Substitute (2.4.22) for FO NCO M224 a7 lL eeanomadds)| ato 


+1 
the subscripts to obtain 


‘ 
Yntjt+1> Me easy - (u — 2d) vr On nee 


(2645259 


Use the coefficients ns to form the linear combination 


a c 


oa fees 


q 
jf 3 ¥nt57 ny og tns jee BEnes-aetl (inal 


a ee hvt. ah (2.4.24) 


0 +) 


Recall tne detinitions (2.4.77) and (2.4.38) cL (Cea and 
Ha respectively and write (2.4.24) in the form (2.4.9) 


to complete the proof. 


Proof of Lemma 2.4.2 


Substitute equation (2.4.2b) in the form (2.4.22) 


aaa EA into, equation (2.4.1c) tojobtain 


rly 


ad ave) ee Kier (2a4225) 
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Equation (2.4.25) is a difference equation in standard 


£ormawhich can be solved to obtain (2.4.10). 
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2 oa bruncatzvon “Error 


The first component Lo Ote 2 1S Used Co minimize 
the truncation error. 
The local truncation error TH of the method (2.4.1) 


is defined by 


a 7 ¥ Tek ac 
ay 8 ST > Gu, TAR Soe eae A ciety 


Yn+1 0 il 


where ve DSethesvector calculated!) from (2.4210) by Assum-— 
ing that past f values are exact and that fc = f(t,y,)- 
The following lemmas express Th in terms of the 


true solution. 


Lemma 2.5.1 


The vector ve Of equation (2.5.1) 1s givensby 


allcaiar eer ae ake! 4? (2.5.2) 
i=1 ' 
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Lemma 2.5.2 


Given swrelation (295.2)7,, equation. (2.54.1) =can be 


written 
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where 


I 
k- 
i 
2 


tal 0, ) 74. (2570) 


Since the error in Yn is at least o(ndt+) each 


+1 


65 an )(2.514)" must equailtione for) 1 “ig «)sSetting 


Xo = Carr = Ul oat (Fall d See U2. D6) 
produces Peel =.) thus, the, local truncatronsergror sin 


the first component of the corrector (2.4.1) becomes 
Oheee ye 


The proofs of the lemmas complete this section. 


Proof of Lemma 2.5.1 


Obtain ve from (2.4.10) by using ve? LnmeplLacemoL 


V£ evaluated at the points t,, ti_jr ---) Cnode2)- That 
is, 
gq-2 : 
Vara etn yma VcmnaA) ees (e507) 
1 aD a = 


Expand each a Ne (2.5.4/) eiedelayLor ser eceobout the 


point t, to obtain 
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canbe bails al ee a uy 
sali hiked 


(a.2.8) . - este - 1) (pyp0"# - 1) = o* 7 


At xozxe aottsgauss LADot oul yeudT .f = tp? esouborg 
asmoved ({.5.8) xodpexi09 ort 20 tnsnogmoo text edt 


ie ae 


-~ 


.aoitose eidt stsiqmon asmmol of¢ to atooiq sAT 
“a , 


to sosiq ni "3eniex yd (OL.b.S) mot Pv mbsddo 
tad? ee Vag? Ap? Bombod sat ts hedeuleve 2 \ 


(v.28), d A aay vi i = 7v an 


edt suods estise 10lysT 5 ai (T -@.8) mi pant’ dos® basqxd 


fiside ot 3 tnteq 


28 


q-2 | 
Tere. tT T 5 
ge Sieo a (fae Leen Gok a 
j=0 
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jJ=0 i=2 y 
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Usinge tnesderintcions (225.3) .06 Oo; reduces (2.5.8) to 
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Proof of Lemma 2.5.2 


Expand ie ee in a Taylor series about the point 


th to produce 
[oe] 
n+1 & ‘ xn, i 


Substitute (2.5.2) and (2.5.9) for v. and hVf_,, in 


2% 5taL)2.to) produces: (285.4) . 
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STEPSIZE AND ORDER 


Se eel Ne rocuce1 on 


In this chapter we outline the technique for 
changing the stepsize and order, examine the error 
introduced by stepchanges, and outline an algorithm 
to be applied when the stepsize is reduced. 

Suppose that at the point th the stepsize h is 
to be changed to h*= ph. Then, the value of the i-th 
component of Xx, must be transformed from hevem yi! CO 
nee sul In order to distinguish between these two 
quantities we denote the latter by xa i° Consequently, 


, 


ie a x* WP fetotel ews ee, Ravel A 
nh n 


* 
we write x) ="colL(x x Bo han 
—1 0) ne; ghd n,qg 


The transformation required to carry x into 
ss is easily determined by using the vector Zz defined 
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where 
= * = 
BO) Seb i De eth), 
which is equivalent to 
Eas 


S22 Stepsize-Changing Errors 


Empirical evidence shows that applying the trans- 
formation R(p) defined by (3.1.2) to x, to change the 
stepsize is a good way to increase the stepsize; but 
substantial errors are introduced when R(p) is used to 
decrease the stepsize. These errors are introduced 
because of the way the interpolation polynomial is used. 
Before examining the interpolation polynomial in general 
consider Example 3.2.1. 

the last, column of Tableq¢3s2.1,shows the error 
as the difference between el ees and Bote) LOG 


=O liye peo, WLC T Cae, (t) is the polynomial defined by 


depp 


x, at¢theppeintjes— b.loggihesejerrorsezare smalljior 


1 
i=0O and i=1, however, they grow rapidly for i larger 
than one. These large errors in the higher derivatives 


of the approximate solution indicate that the solution 


polynomial oscillates about the true solution. 
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Example Sey aa 


Integrate the initial value problem 


y = Cy iy lhl, Ca eB) 


one step forward’ with h = 0.1, gq = 3, % = ¢ol(3/3, 1, 
3/45, 270). tCalculating the derivatives of fthe’ true 


solution y(t) = t? yields the starting vector 


Xp = COMCE Ost Uie pads. OL) (3 e204) 


which defines the starting polynomial 
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0.01) defines the polynomial 
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column Of Table 2..3.1, defines the polynomial 


P, 1 (t) = 1.6105125 + 7.3205(t-1.1) + S39 (t-1.1)° 
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Values ofhthe Polynomiale+ 3.253 )itand ,(3.2.4) 
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We shall now examine how the interpolation poly- 
nomial in Nordsieck's method is used when the stepsize 
is changed, We derive a transformation S which takes 
the vector x, EnLOeanOtneresvector Pa which represents the 
interpolation polynomial. P(t) in terms of past values. 
The transformation S and the vector p,, are USeCascLOetLans= 
form Nordsieck's formula into an equivalent formula using 
past values of P(t). We then show how the oscillations 
present in P(t) are translated into the scaled derivatives. 
It can then be seen how step-changes introduce errors 
into the solution vector x 

Recall that SG Consists ol thes 7—ch scaled 
derivative of P(t). Expand Es co Ney ahvee fe) URE AR ep 


series about the point t, to obtain 
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1 0 0 0 0 
0 1 0 0 0 
S = }|0 elem 3 2 3/74 — 17 : 
0 WF pes, 570 27 aal/aG 
0 1/24 7.8 1/8 yee 


To find the equivalent of Nordsieck's formula 
in terms of past values write equation (2.4.1b) in the 


form 
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equation (2.3.5) to produce 
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(i) No starting errors are committed. 

(ii) A constant stepsize is used. 
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Suppose now that at the point ce the stepsize 
h is increased to 2h. Retaining conditions (i) and 


(iii) above, Be is then replaced by xe where 


R(2) 0x 
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Relation (3.2.10) shows that Nordsieck's method does not 


use values of f to continue the integration, rather, 
it uses values calculated from the interpolating poly- 


nomial P(t) at points spaced a distance 2h apart. 


These points are marked by (t) in Graph 3.2.1. Examin- 


ing Graph 3.2.1 more closely indicates that the errors 


_ ee) ae 
ones foi P (Ces are zero for i=0,2,4 and are 


non-zero and of the same sign for i=6,8. Moreover, 


these errors €_ . increase in magnitude as i increased. 
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However, premultiplying Py byethesmatrix |Setoslorm the 
vector te causes, due to the alternating signs of the 
row entries of S, these errors partly to cancel. In 


addition, the row entries She strongly decrease in 
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Graph of P(t) and £ versus t 
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magnitude with increasing j when j > i. This reduces 
the contribution of the large errors. Consequently, 
no difficulties arise when the stepsize is doubled. 
Clearly, a similar illustration and conclusion can 
be given for increasing the stepsize to an arbitrary 
h* using a method of order q. 

Suppose now that the stepsize is halved from h 


tOen A2aeitnenin 63.«2,5.80 isbecomes 


35a us P. (32 ol) 
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In this case values of the interpolating polynomial 
P(t) spaced a distance h/2 apart are used. These 
Values face\macked oc) eineGraph «3.221. Ut, can nowpbe 
seen from the graph that some of the errors which are 
close to the point t, are non-zero. These errors are 
then jmuLtipliedsebysatheslarge.cpiniess#oh.S.. .bhus,~rela— 
tively large errors in xh result. Consequently, care 


must be exercised when reducing the stepsize. 


3°3, Tests for Stepsize and Order 


We recall that in section 2.5 the parameter Loy 


is determined so that the (qt+l)-value multivalue method 
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ane Unfortunately, 


has a local. truncation error O(h 
2iathegerror as 0 (not) then good estimates of the 
Locals truncation errors: are difficult to obtain. 
However, if the error is o (natty then good estimates 
are obtainable. 

After experimenting, the following method of 
choosing 2% 1s used. Let Ea be the vector which makes 
the multivalue method (2.3.1) equivalent to the gq-th 


order Adams-Moulton method. Our program uses vectors 


2 defined by 
a9 


we 
HT} 


ard Lokw iencds >o0 


This choice insures good stability characteristics and 
produces local truncation errors O(n 1a 

To see why it is desirable not to use maximum 
order (q+tl)-value methods suppose that such a method 
1s used? Thestruncation error of the method is out oe 


Then, truncation error estimates for the (q-1),q, and 


(q+1l)- value methods are required. These estimates are 


(qt+1) (q+2),qt2 
of the form eae Mees) “ee yet h ; 


(q+3) mots 
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h ae and 


respectively, where C G 


qt Seats 
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However, in practice we found that os Fel approximates 
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Consequently, backward differences of de ee are not 


/q! with only one significant figure accuracy. 


accurate with any signifcant figures. Therefore, two 
of the required estimates are not accurate. 
On the other hand, if the q-value method has 


truncation error Ont} then the truncation errors 
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where C an Chec are known constants. In 


Gu Cot1! ehiaa 


practice, we found that ries and a Fel produce 
estimates which are accurate to at least one signi- 
ficant figure. Thus, with this choice of the parameter 
Ro two of the required three estimates are reasonably 
accurate. 

USing the. estimates’ (3.33.2), (3.3.3), and (3.3.4), 
of local truncation errors we calculate three stepsize 
estimates hoew 


using these stepsizes and formulae of corresponding order 


hoe and h* respectively, so that 


qriy 
produces local truncation error estimates which are equal 
to a given absolute local truncation error tolerance 

ie 
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Lovecalculate ho-l wemuses the fact. that the wlocal 


truncation error for the method of order (q-1). has the 
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The squantucy, hoy is now an estimate of the stepsize 

which would have produced a truncation error estimate 

equal to ise had the formula of order (q-1) been used. 
In exactly the saméejway “353.3) and (°.8.4)* are 


used to calculate the remaining estimates 
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| In*practice /'‘these*choices’ tend'*to cause too 
much stepsize changing. Thus, when changing stepsize 
the estimated stepsize should be reduced by a factor 
a < 1. Whe »pprogram.jin, the, Appendix, uses. a. = 0.9, 

We must now decide which De to use, The algo- 
rithm which follows shows how the program in the 
Appendix chooses the stepsize and order. 


(35) Specify a set of bias factors Bg-1' B +1! 


and 
g Bg 


which must be chosen so as to reflect the amount 
of work expended when using order q-l, q, qtl, 
formulae respectively. 
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° * * ; 
(iv) Calculate eat Set and Pot] ees 7 cand the 


number q* so that 


* 2 * * * 
ioc io Bluth: char candiens iumegairecsle axed Ra 
set h~= one hoe peande cine Ordere longa. 

(v) Teh 23h do, (Vii) next. Lf therstepsize is 

reduced four times on a given step reduce the 


Orders to.2: 


(vi) Repeat gcneysctepewitn. SCepSi1 zen andene turn eto 
Cals) a 
(vii) Set the stepsize to h*. 


(viii) Accept the step and return to (ii). 


sya) De Stepsize Changing Algorithm 


In this section, following a suggestion by Krough 


[14], we develop the correction 
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Note that before changing the stepsize to h* a 


tentative step is made with stepsize h, therefore, 


Goat is*available without ‘extra function evaluations. 
In our analysis we assume that 

(eis) The exact solution has been obtained at the 
POIncs tos tie SE A She it 

(i) The solution is a (q+tl) degree polynomial. 


(iii) A gq-th order Nordsieck metnod is used with one 
application of the corrector formula. 


In addition, the predictor matrix will have dimension 
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Furthermore, we augment { With, a zero,.so that, 


L= col(£o, Qarerer & oy 


gq! 


Making two steps with Nordsieck's method yields 
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Ene enmasy Of “section 3.2.5 Using equation (3.4.5) and 
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CHAPTER LY. 


OUTLINE OF THE COMPUTER PROGRAM 


ALIS Uintroduction 


A computer program, called NDSK, has been written 


to solve the system of ordinary differential equations 


(1 
Bh ee eee y ey); ep ee (44-1) 


y; (to) PM 


where 1 < Is 10, andt, <s ts ty. 

NDSK uses Nordsieck's method to generate a sequence 
of solution points ¥(tj), ¥(ts), ey Y¥(t,) on the basis 
of parameters specified by the user. At each step the 
stepsize and the order are determined by the program; but, 
provisions are made to enable the user to override the 
programs choice of the stepsize and to limit the order. 
This ability to vary the order allows NDSK to handle mild 
discontinuities in the driving function £ and its deriva- 
tives. 

A general description of NDSK is provided in 
section 4.2 and details are given by comment statements 
ineither program Listing. * Sections 4.20 provides! thetdetaiils 
of NDSK's calling sequence, and section 4.4 gives the 


calling sequence of the subroutine FEVAL which evaluates 
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the driving function. Section 4.5 provides some detail 
on how the stepsize and the order are chosen. 

Since the present version of the program is the 
product of considerable experimentation, where we have 
tried several innovations that appeared promising, we 
conclude this chapter with a brief outline of its deve- 


lopment. 


Se 2eucneral Description Oma~tie Program 


The user's calling program communicates with 
SUBROUTINE NDSK which organizes and supervises the in- 
tegration process by calling several subroutines that 
perform various functions such as: 

(i) calculating the next solution point, 

(12) testing for stepsize and order, 

(iii) defining constants and coefficients, 

(iv) printing the solution. 

The interaction of these functions is indicated in 
Figure 4.2.1. 

When control passes to NDSK the input arguments 
are checked for errors. If any errors are discovered 
these are corrected, if possible. The arguments are 
then transferred into the working array A(10,18) located 


in COMMON BLK1 which is available to all subroutines. 
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Next, NDSK calls SUBROUTINE START which initia- 
lizes another part of A by calculating starting values 
required by Euler's method. Also, START initializes 
two vectors, CP(12) and D(10), to contain the error cons- 
Lante tore thelg-theordéerimethod(iniCPla)yifandito Contain 
the coefficients ge in D. Furthermore, START is called 
whenever the stepsize tests fail four times on one step. 
Such repeated failures indicate an excessive error build- 
up in the higher derivatives. In this case the higher 
derivatives are discarted and the integration is started 
over with Euler's method. 

After the initialization, NDSK calls SUBROUTINE 
NORD which saves the values of A in another array ASAVE 
and advances the integration one step forward with step- 
size HN. NORD is also used to interpolate the solution 
tosthenprantspeintsoandptne endpoint. 

Soon as control passes from NORD, NDSK calls 
SUBROUTINE TEST. When entered, TEST calculates a new 
stepsize HW and a corresponding order on the basis of 
thenalgorithmrgiven, inAsectionys. 820, LETHWOOHNKDtEheCstep 
will be repeated with, perhaps, lower order; otherwise, 
The step 1s accepted . 

The sequence CALL NORD - CALL TEST is repeated 
Hinclieonelcr. che, LOLLOWINngsconcd1 tions 5 salisraed ; 

Cas) The stepsize HN is reduced to a value less than 


HMIN. 
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Cin) A print-point is reached. 
(111i) The endpoint is reached. 
At this point NDSK calls SUBROUTINE PRNT to initiate 
printing the solution value. In case (i) PRNT prints 
a message. In cases (ii) and (iii) PRNT calls NORD, if 
necessary, to interpolate the solution to the print- 
point. When control returns with condition (ii) NDSK 
returns to execute the sequence CALL NORD - CALL TEST, 
otherwise, NDSK returns control to the user's calling 
program. 

Finally, we collected all parameters in SUBROUTINE 
COEFTS which is called by NORD whenever the order changes 
to place the corresponding parameters in the arrays C 


and D. 


4.3 The Arguments of NDSK 
To integrate the system (4.1.1) of differential 


equations the user must 
CALL NDSK(TO,TF,TPO, HMAX, HMIN, YO,NEQ, ECT, IODR,KDEL,CQ), 


where the parameters have the following significance: 


TO. is a, real variable equal to the initial, value.of 


the independent variable. 
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NEQ 


ECT 


mE 


is a real variable equal to the final value of 


the independent variable. TF<TO is permissible. 


is a real variable equal to the spacing between 
print-points. The solution will be interpolated 


anos pulicedmatetnespOolnesmia =" TO me Ke lPO meK=0 lee 


k 


where L depends on the value of TF. 


is a real variable used as an upper bound for the 


stepsize. 


is a real variable used as a lower bound for the 
stepsize. HMIN should be small. In most cases 
its value should not exceed 10°2°, for, the 
program starts the integration with Euler's method 


which must use a very small stepsize to produce an 


accurate solution. 


is a vector of real variables containing the ini- 
tial values of the dependent variable. YO should 


be dimensioned to 10. 


is an integer variable specifying the number of 


equations in the system (4.1.1). 


is a real variable used as error tolerance. ECT 
appears also in the parameter list of the driving 
program FEVAL, hence, its value can be changed at 


each step. 
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is an integer variable specifying the maximum 

order to, be jised. 

is an integer variable specifying the minimum 

number of steps that the program is to execute 
between order changes. 

is a vector of three real variables related to 


the bias factors 8 eis Tnasecuiones«o. 


Gat ee qa. 

That is, we use CQ(I) = By. I=q-l,q,qtl. These 
values may be used to bias the program towards 
using lower order by assigning a larger value 
bonCo (5) stiang £6. CO 1) vand CO (2) 2 vAwreversal 

of the relation between these factors introduces 
bias towards higher order. Normally, values in 


the interval [1.0, 1.5] are used. In most cases 


COPS te pls eo WOLKS wells 


The variables; TO, TF, TPO, HMAX, HMIN, and YO, shouid 


be DOUBLE PRECISION with YO dimenioned to 10. 


4.4 


Dngeto 


Thee Oey Lig mune clone liVAT 


The user must provide a driving program conform- 


SUBROUTINE FEVAL (T,Y,F,TO,TF,HN,ECT) 


where the parameters have the following meaning: 


is a real variable containing the current value 


of the independent variable. 
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x is a vector of real variables containing the 
the current values of the components of the 


dependent variable. 


F is a vector of real variables into which the 


value of £ must be placed. 


TO,TF are as in NDSK. The ability .to change these 
values is useful to restart the integration at 


singular points. 


HN is a real variable containing the value of the 
stepsize used to make the current step. Speci- 
fying HN overrides the programs choice of the 


stepsize. 


ECT is as in NDSK. For many problems it is useful to 
specify an error tolerance in terms of Y or F. 
Several examples of such practice are given in 


Chapter sv? 
Thes varlables TY F710, 1,F;, and HN, must. be 


DOUBLE PRECISION with Y and F dimensioned to 10. 


4.5 The Tests 


The stepsize and the order are determined by 


SUBROUTINE TEST . 
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To calculate the stepsize estimates, TEST expects 


the highest order scaled derivative of the I-th compon- 


ent of the system (4.1.1) in A(IOR). Test uses the 
variable IOR for q which is used in the thesis. Fs 
and VE 41 are expected to be in A(I,14) and A(I,16) res- 


pectively. Moreover, the square of the inverse of 
(t/q): uscd in=equation (323.5), and the factor gq? muse 
be in CP(IOR) and C(2,10R-1) respectively. Then, TEST 


calculates the factors 


1 
NEQ 2 2*10R 
T1=CQ(1)*(CP (IOR-1) #c7(3,IOR-1)* J Cae ; 
ca) 
D 
NEQ 2* (IOR+1) 
T2=CQ(2(CP(IOR) *C*(3,T0R-1)* ) (ATA?) 
T= 
i 
NEQ 2* (IOR+2) 
T3=CQ(3)(CP (IOR-1) C7 (3,TOR-1)* J (ALTO) *) 
I=1 


where WTS(I) is the maximum value of the I-th component 
of the approximate solution for to Sey Ss to: TEST uses 


these factors to calculate the new stepsize 


HN 


HLD/RHO , 


where RHO MIN (T1,T2,T3) and HLD is the stepsize which 
has been used. The order is then adjusted using the 


following rule: 
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(1) If RHO = Tl and IOR is greater than one the order 
is decreased to IOR-1l. 

Cy) If RHO = T2 then the order remains IOR. 

(iii) If RHO = T3 and IOR is less than IODR the order 
is increased to IOR+1, where IODR must be given 


in the calling sequence of NDSK. 


4.6 The History of the Program 


The first version of our program used the method 
as given by Nordsieck [16]. This version restricted 
stepsize changes to the form h* = ooh where k is an 
integer in the range -5 to 5. It used fixed order methods 
Loecicerandes | slog eanGerteusecastierstalting procedure 
outlined by Nordsieck. This program used a very conser- 
vative, often extremely small, stepsize. 

We then implemented facilities to allow stepsize 
changes of the form h*= ph where 0.1 < p < 10.0. At the 
same time we changed the starting procedure. The new 
procedure started the integration with Euler's method 
and increased the order on subsequent steps. These two 
changes improved the efficiency of the program consi- 
derably. 

Next, we implemented facilities to let the program 
change the order. Initially, we had the order-changing 


mechanism biased towards using low order and used the 
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Vector si. as given by Nordsieck y[16]/) This#version 
experienced two difficulties. Firstly, the local 
truncation error estimates were poor and, secondly, 

the order was changed very frequently. Actually, these 
two difficulties are related. 

| The final version of our program allows the order 
EORV aL moe tween 2eandel2 elt uses the parame cers (las 
given in section 3.3 and uses the algorithm given in the 
same section to change the stepsize and order. This 
algorithm biases the program mildly against order 
changes. The algorithm of section 3.4 is also imple- 


mented. 
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CHAPTER V 


RESULTS OF THE TEST PROBLEMS 


5Spaci Introduction iand Summary 


The problems are intended to cover a realistically 
wide spectrum of problem type. We included some fairly 
difficult ones, some with singularities in the solution 
or in its derivatives, some with oscillating solutions, 
and a class of stiff systems. Most of these problems 
involve functions which are relatively easy to evaluate 
and whose analytic solution is available. 

When using our program to solve a given differen- 
tial equation an error tolerance must be specified. In 
this context it should be noted that we restrict the 
absolute error rather than a relative error. This 
choice avoids difficulties when the solution is near 
zero. Moreover, we restrict the Euclidian norm of a 
vector of local residuals. Alternatively, facilities 
can be implemented to assign error tolerances to the 
individual components of systems of equations. Further- 
more, provisions are made to allow the user to specify 
the error tolerance in terms of the solution of a given 
equation or in terms of its derivatives. 
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guidelines can be given. In general, the error tolerance 
should be small whenever the error grows strongly and 

it should be relaxed when the error grows weakly. More 
Specific, if assolution accuratetto gq significant figures 
is required then an error tolerance proportional to 

107 (att) will, in most cases, produce a good result. 
Furthermore, the error tolerance must not be smaller 
tiene chesulnicertainty n= tne urltving = tunectron. =what=1s, 
if the value of the driving function is given accurately 
to a decimal places then the error tolerance must be 
larger than 10°”. 

For particular examples of choices of error tole- 
rances we refer the reader to the examples. 

Comparing results obtained with different programs 
is, in general, not easy. In particular, such comparison 
requires a basis which is difficult to define. However, 
we believe that the number of function evaluations 
required to produce a solution of specified accuracy 
AS agvalidscriterias),On,this basis»wep compared our 
results with those obtained by Gabel [4]. We found that 
our results were obtained slightly more efficiently for 
equations having solutions which decrease in absolute 
value, but the reverse was true for equations having 
increasing solutions. Moreover, we observed that our 


program increased the order and stepsize slower which 


results in a slight overall Loss of etiiciency. 
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Bi Equations of Exponential Functions 


The two simple initial value problems 
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which have the analytic solutions y = a and y = ae 


appear quite often in reports on differential equation 
solvers. 
We obtained solutions at t = 20 using, for 


equation (5.2.1), the error tolerance 


ECTs 10am (Som) 


and,) Lor cquation. (5.2.2), 


ECT = 10 eK lo sis Pee) a (5.2.4) 


it 


The results are shown in Tables 5.2.1 and 5.2.2. 
Equation (5.2.2) is an example of the usefulness 

of variable error tolerances. To see this, let us 

assume that (5.2.2) is solved using an absolute error 

tolerance ECT = Ware, Such bound on the absolute error 

together with the initial condition y(0) = 1.0 should 

produce, near t = 0, a numerical solution with about 5 

-8 


Significant figures accuracy. However, y(20)) = 10 9, 


hence, using the same absolute bound near t = 20 cannot 
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be expected to produce a solution at’ t ="20"having any 
accuracy at all. Therefore, we use the error tolerance 
(5.2.4) which is proportional to the magnitude of the 


solution. 


5.3 Sine-Cosine Differential Equations 


The initial value problem defined by the system 


of differential equations 
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an upper bound HMAX = 5.0, and a lower bound HMIN 
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table 5.3... shows the absolute errors aa 4 
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i=0,1 and the number of steps taken from 


eae 
12 


= 0. For simplicity we show ee ; * 10 ; 
v 


Results of the Sine-Cosine Problem 


STEPS 


10 OS As7i83 Oe 2he0 81 
218) =O O90. = Oh 20 282 
100 = ier? O39 Qezou7 534 


Oro Zon = ie Ou) A 


DATS bie. Oren 


The program evaluated £ about twice for each 
step taken, namely, 2106 times over the interval 
0 < t <¢< 200. The relatively large number of steps 
taken to reach the point t = 10 indicates that the 


program increases the stepsize and order slowly. 


5.4 A Non-Linear Problem 


The initial value problem defined by 


=A RY, F St thd 


63 


ny 


a . Per otk 


= 4. ,2 etoxms etuloeds. oitd awore f. €.¢ po 


mex? asxe+ eqote to sedmun oft bas f0ek tiegeee™ tat 
"OL |, 3 wore sw yoiorlqmie 10% .O=2 
_ 


moldoxwd oentead-snLe ads to esivasd 


L.t.¢ Saar 


Hoss xot soiwd tuods 1 bedsulsvas maipoxg SAT 
fsevissani sds xsvo semis QOS Di ald ged 

ageje to xusdmwa spzsi ylevissiex sont 00S 24. a) 

end gedt astsotbat Of = 3 stnioqg erit ddpex od igdad 


vy” 


-Yiwola xsbio Bris Ssteqsta sz ssrastoni msirperq 


nodoat sag —— 


a 


a) : | sitll a 
; - aeon: asisoug aulev fa 3 : od? 


7 - a . 7 a 7 en 


-_ 


“F ‘au - : | 


64 


> 
II 


VC) e = coll 070. 0,750. 077 1,0)n; 
is a non-linear system which has the analytic solution 
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respectively the radius, the radial velocity, the azimu- 
thal angle, and the azimuthal velocity of an orbit. The 
units are feet, feet per second, radians, and radians per 
second. 
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The third component should attain 


Weasolwved problem (575.1) over 10” orbits using 
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y(0) = 1, and 
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y(-4) = 4. 
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reduced the order to that of Euler's method and the 
stepsize to aarp before passing the singular point. 
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The results are reported in Table 5.10.1. ‘We 
noted that a relatively large number of steps were taken 
to reach the point t = 10 which indicates that the pro- 


gram had difficulties starting the integration. 
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We note, that«the derivatives of.equation (5;11.1), 
starting with the third, have a pole at t = 1. There- 
fore, the order of the integration formulae must be 
lowered in order to pass the point t = 1. Equation 
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We obtained solutions for equations (5.11.1) to 
(5.11.3) 20versthe interval 0... jt<.2 anddfornt(>elhls4):tover 
the interval - 5 st < 1. We used the following error 


tolerances fot .equations.a(5 sl lnk) sto, (Golle4 peespectively: 


Bene Shy it Pa (5.15) 
-10 ary eee 
ECT = 10 REM ERETG c/a, (or. eso) 
n n 
Senet “sik ae 
BC Olas) cea| ete Lo ae, (IBS a) 
ECT = 10° . (5.11.8) 


The results shown in Table 5.11.1 include the numerical 
solution at the endpoint, the error relative to the true 


solution, and the number of function evaluations NOF. 
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ResULCSsO ime ual tons oil. 1) ttOul oem. 4) 


EQUATION Y(t) RELATIVE ERROR 


Gory a) 


Ost 23057126 O72 .33 /40-06 


0774999959 


(Ou ae) 0.13: 39ER—-06 


(Cone lees) Des Oe oes: -0.1441E-03 


(ot 134) US LG Ag/ 


0.83/71#H-08 


TABLE S24... 


5.12 The Hull and Cremer Problems 


Hull and Cremer [10] list 17 equations which 
they used to test cine efficiency of iprediceer-Ccorrector 
schemes. These equations are listed in Table 5.12.1. 

We solved these equations and report the results 
iiptablewo.i2.2.  elneeanitials condLGionsmet t= 0 were 
obtained from the analytic solution. For each problem 
we list the error tolerance ECT that had been used, the 
upperbound HMAX for the stepsize, the error relative to 
the true/solution at t = 40, and the number of function 


evaluations NOF. 
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The Hull and Cremer Problems 


PROBLEM 
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Results of the Hull and Cremer Problems 


RELATIVE 


053162E-05x | Y| 0.1950E-04 
0.1414E-04x|yY| 0.5958E-06 
0.1414E-04x|yY| “0.8912E 14 
0.3612E-04x|yY| -0.7510E-06 
0.2718E-04x|¥Y| 0.7304E-04 
0.1E-04x|¥Y| -0.4828E-01 
0.1E-04x|yY| Oo s ones 
0.1E-04x|yY| -0.1832E-03 
0.1E-04x | Y| 0.1128E 01 
0.4E-05x|¥| —Oseecln-04 


0.1414E-04x|yY| 0.1084E-05 


0.JE-04x|yY| 022137203 


0.1E-04x|yY| -0.1130E 15 
0.1E-04x|Y| -0.4624E-05 
0.1E-07x|¥| -0.1019E-04 
0.1E-04x|¥Y| 0.1034E-07 


0.1E-05x|y| -0.5024E-05 


TABS 45. 12. 2 
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Four of the equations listed above, these are 
starred, are unstable. Our program integrated all four 
of these problems to the endpoint. 

Our program and Gabel's were run on different 
machines, hence, comparisons. of.results are difficult 
to make. Yet, we have observed consistently that 
Gabel's program produced better solutions for increasing 
functions, whereas, ours produced better results for 


functions decreasing in absolute value. 


peo Seen oOLUCLOnNEOn y?? = y log y/t 


The family of curves y = aE Satisfying scne dit— 


ferential equation 


y‘)) nV LOGY /.t (Saas el) 


nas a icommon gpoimtygat th= 0.72An initial condzbion, given 
at any point t # 0, uniquely defines a particular curve. 
However, due to the common point at t = 0 and the exponen- 
tialecormect thestami ly sOresOlLuLiOn Curves errors =1ntro— 
duced, nédwete—-0sare.greatlyomagnitredss Hence, near the 
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Note that ECT attains a minimum when ae ieeenat 174 


when t = 0. The results are given in Table 5.13.1. 


(1) 


Results of y Say MLOC ayy tC 
ie RELATIVE ERROR STEPS FUNCTION EVALUATIONS 
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det4e Lie oOlutionlol a Bessel equation 


To test a differential equation solver over a 


long interval Nordsieck [16] used the equation 


£29 (2henl ty (huapade owes?) yi th \Gak (5°45 45 


where y(t) = J, 6 (t) which is a Bessel function of order 
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loz Ji (t) is small in absolute value at t = 6. How- 
ever, its value increases rapidly when passing t = 6 
and oscillates more than 1000 times over the interval 
6.55 Sboleo: 


We integrated (5.14.1) in the form 


PU RENE aol vaio) = 0 1201950% 1048.3 
256 ~ 
oe (=3- Vy -¥y/t , y,(6) = 0.2986480 x10” , 


using the error tolerance 


HGT @= LOG! ({up¥A || + heda? fap) 


We specified HMAX = 5 and HMIN = I = UltS Geeisulias Eig 


SHOWN LAD lem je aie 


Results of a Bessel Equation 


: FUNCTION 
| | __srron jones PUREE! 
-0.009745831|-0.000001150} 23869) 69464 


TABLE 5.14.1 


6136 |-0.0097446808 


The results compare well with Gabel's [4]; but are 
not as good as Nordsieck's [16]. Actually, Nordsieck 


obtained his solution using stepsizes of the form ae where 
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were 2928 > 


mzoi off ni (LA1.2) botexpesik) ow 


Rapes 


POL x d2eross.0 = (2) py s° gta 


oy - . - i) VP (f) 
, FOL x OanaBes.0'= (Ol py « N\gk- o¥tl = = Oy 


sonsielod totze end pr. 


i Ge Te 
aan || + || Mp jb OL = foe = 


3 ' 
Sun atives: sat vor 2 WIMH Bos @ = %AMH beiiisegs sh 


.f. 81.2 olds? at awore | 
? 


dorseupd Iseesa s ae etivesh 


=o 


L.sf.¢ ataay 


x8 gud <1) ea'isdso div flow STAgmoD. etivesx oat 
Yeeiebxow \yiisuas4 .|21) 2'tostebuew es boop as 
| ee eA nd ad ‘poten nokvvLo8 ait & 


81 


K is integral. Such choice of stepsizes reduces the 


round-off error considerably. 


5.15 The Solution of y"?) = a*b/(t*+a’) 


Selecting HMAX for problems of the form 
NOD are aS) pb 


ys "quite important; partrcularily; "if ais’ smali* Equation 


(5.15.1) has a general solution of the form 
Vite anarctans(t/bD) rt a arctan (-t)/b) 


Witch “as ea-qrapleag-r lustratedsin Figure: 5.15.1. 
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This example involves searching the t axis for an 
extremely narrow interval in which £ peaks. If HMAX 
is larger than this interval then the program will 
bikelyyoverstep it without noticing the sudden change 


and produce a useless result. 


We solved this problem with a = pt ese, eae oe 
for - 5 Sa SS =. We used y(-%) = 0, error tolerance 
me S yy "| rel EE S Oo egitevers ieee) alah hil 


The results for each value of K together with 
the solution Y(%), the error at t = - the number of 
function evaluations NOF, and the number of steps, are 
Given ere beaolew oo. Lr. 


The analytic solution is y(%) = 0.37450728 x TO e 


Results of Equation (5.15.1) 
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5986 “The (Solutwon SoffODE's “with Eroresstimates 


Gabel (4) gives an crrom analysis .Or the initial. 


value problem defined by the simple system of differential 


equations 
a Feel ees 
ii 
y: gs AGS 5) 
1Y5_4 Hie ee La) 79.7 oF cee ete ay 
Yah} = 0.0 Ihessysctem. (5.16.1) has the general solution 
y(t) = t*. The analysis concludes that if an absolute 


error tolerance t is used then the error oan F in the 
lA 
humerical solution,for the i-th component of .(5.16.1) 


should obey 


(Sy iSong Rae eye cinh ys (Salon 
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Moreover, if a relative error tolerance t, . = T|Y_ .| 
n,i n,i 
is used then 
Sa ue ele ye 
epee Sake yAlG awe) (se lGes) 


ought to hold: 

Werintegrated (5.16.1) using K-86, Eno wan and 
allowed the order to vary up to and inclusive the value 
9. Then, we repeated the calculations with the order 


restricted to values less than or equal to 6 which is 


less than the degree of the solution polynomials of the 
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tases twOorcomponenes Of *(>216"1) > FTabtes*™5o16r1l4ane 5.46.2 
list for i=2,4,8 the value of the analytic solutions Y, at 
f=? andst=20)} the value ot the bounds (5.16.2) and (5.16. 3) 
using the final value of t, the absolute error, and the 
number of function evaluations NOF. 

The results in Tables 5.16.1 and 5.16.2 show a 
much smaller error than bound where iis small. However, 
when 1 equals 8 the errors are approaching the values of 
the bounds. This phenomenon is due to the fact that our 
program does not bound the errors of the individual 
components of a system of equations; it bounds the 
Euclidian norm of the vector of residuals. Consequently, 
some variation of the error over the components of a 


system of equations is possible. 


RESULCSEOL -EQUACIONS (5...0.1)), K=—.o, ORDER < 99) at, & = 12 
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RELATIVE TOLERANCE ABSOLUTE TOLERANCE 
y (20) 
ERROR BOUND ERROR BOUND 
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5.1l/ The Parabolrvegsurtace 


An example, given by Gabel [4], to test a diffe- 
rential equation solver's ability to handle discontinui- 


ties is the system 
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Note that the parabloid vs a" a = t separates the 


Euclidean 3-space into two regions Ey and E, as indicated 


in rigure 5.1/.i and that any solution of the, system 
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(o.27.1)Ssbharted=in Ej, except Yor ye 0, takes a 


finite jump when meeting the surface ve 25 vt = t which 


from E i 


separates ED 3 


PILGURE Oey ail 


We attempted to calculate several solutions of 
(Sieol/-L)pusing-—varlous anitial values at t = <1... Our 


routine had no difficulty extending solutions across the 


mentioned surface whenever ve + i < >. However, all 
solutions were quickly terminated whenever ve +yt 2 >. 


Some of the results obtained are as follows: 
(i) Using the initial condition Yo (-1) = y, (-1) = 0 
leads to the trivial case Yo (t) = y, (t) = 0. 


The calculated solutions followed the line 


L.VLi,e Gavort. 
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(ix) Using Yo (-1) = y, (-1) = 0.25 the solution curve 


enceread gener regron fk wirrout rrr culty; but, 


2 
the curve terminated when it touched the surface 
from inside E, at the point (t,Y¥o,Y,) = (2° S405 729, 
“Oeoua/) Dp Leola ee At Sires poire = Waa been 
evaluated 239 times. 

Cire, Ustirg Yo (-1) = y;(-1) = 0.5 the solution curve 


encountered the surface at (t,Y Y,) =S(0ES 7625 


On 
OT IZ2022 ye OT ool) ana terminated. The Luneceion 


£f had been evaluated 70 times. 


5.18 The Ricatti-Bessel Equation 


The Ricatti equation 
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has the solution 
PU EO Gi, (ee 72) 
3/4 -1/4 J 


where Haye and Beat ye are Bessel Pan Cen Ce Oe orders 
3/4 and -1/4 respectively. 

The function w(t) has a pole at any point t. such 
that Da ca 2) = 0; hence, a differential equation 
solver should not integrate past such points t.. There- 


fore, it can be used to locate the zeros of the function 
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co ee To locate the zeros of both functions, 


ie U pita vel and he ngh abi we used the algorithm: 


(i) 


(11) 


acahaty) 


(iv) 


the error tolerance ECT = 10 


Starting at oe t. integrate the equation 
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( 5nd] Seo XG vs at which point the program 


sic all 


reduces the stepsize to a value less than the 


specified HMIN. 


Take the value z. = ee 


i+l 44172 aS an approximation 


of a zero of Shae i ah ced 3 


Switch to the related equation 
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which is obtained from (5.18.1) by applying the 
transformation u(t) = w(t). 
Staygtingvater = t. integrate equation (5.18.2) 


tO Ghee point c= 9c. where the stepsize is 


1+2 


again reduced to a value less than HMIN. Take 


Zz = Cee to; be an approximation of a zero 


L+Z 


ore ord] Replace t, by tia and return to 


See 2 


step (i). 
Using the above algorithm with HMIN = 10 ~ and 


i we calculated approxima- 
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Tables 5.18.1 and 5.18.2 list these zeros together with 


the true zeros y; which are obtained from the tables ana 
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The large number of function evaluations is due 
to the facts that the integration must be restarted 
whenever a zero is located and that our program has the 
tendency to reduce the order and stepsize before giving 
aos 

Table 5.18.3 shows ye numerical solution w(t) 
at the points) BeewK) Ke b,2770,4,5. whe true values 


w(t) are taken from Gabel [4]. 
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Results of the Ricatti Equation 
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Belo Arr Oscillating Solution 


The initial value problem 


. -2ty, log Sey ne Yq (9) =erevy 


| 
fH 


’ Oar eels) 


ke 
aoe 
I 


2ty, log yg + yj, (0) 


is used by Fehlberg [8] to test Runge-Kutta methods of 
Orders] lLOuGc.) eine system (5.19) Nase the analyte 


solution 


EXP (cos as ; 


Yo (+) 
vy, (oRRSGexp (saint), 


hence, both components of the solution vector oscillate 
with period /n and amplitude Bala. 

Fehlberg [8] uses a local truncation error tole- 
rance equal to ‘jee In our environment such small 
error tolerance introduces serious roundoff error problems. 
Hence, we experimented to determine the smallest values 


which when used as error tolerances allows our routine to 


perform reasonably efficient. The two values 


ECT = 107° (5.19.2) 


and 


ECT = TO aan tee by tom 
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satisfied such requirement. 
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Wepintegrated (5.19.1 )Mover the anterval Oe<'t <5 
and report, the results in. Table 5.19.1. 

The Verrorein Ourasolucions iS. largerethan sin 
Fenlbergeas bysartactor ofvabout 10.) However, wercounted 
6033 and 7186 function evaluations using reupect ais 
the error tolerances” (5719.2). and (5.1923) Gwhichras con— 
Shusieca iy below the values reported in [8]. Thus, this 


result supports the conjecture that multistep and related 


methods are more efficient if £ is expansive to evaluate. 


Results of an Oscillating Problem 
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5.20 An Unstable Problem 


To test our routine with an unstable problem we 


integrated the system 


yi’ =Ayte CEO; (eo 20a.) 
where 
0 ils 0 0 
eee 0 0 uf 0 : 
0 0 0 a 
ih 0 0 0 
We used the initial condition 
Va) eC OL Oper 2 3p 0 pees eee (522022) 


The ;differential system (5.20.1) is equivalent to 


the equation 


which has the general solution 


ea —t it —t 
y(t) = c,te + coe + c3e + c,sin(t) t+c,cos(t), (5.20.3) 


where the c,'S, 1-1,2,2,4,>, are constants which are: de- 
terminedsby the initial condition. | in our case (5.20.3) 
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which is free of terms containing ane However, small 

errors generated by the numerical method introduce such 

terms which eventually dominate the required solution. 
We integrated the system (5.20.1) over the inter- 

Val OV CC < 40 using the error tolerance ECT = Roe 

The results, listed in Table 5.20.1, show the error in 


the first component Yo 7 the BucIi dean norm OL the errors 


over all four components, and the number of steps. 


Results of an Unstable Problem 
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CHAPTER VI 


CONCLUDING REMARKS 


This thesis consists of a study of Nordsieck's 
methods for solving numerically initial value problems 
in differential equations. The method is developed in 
Chapter II. The stepsize changing is examined in 
Chapter III. A brief outline of an implementation of 
Nordsieck's methods is given in Chapter IV and the test 
results are reported in Chapter V. These results support 
the conclusion that Nordsieck's method definitely deserves 
consideration as a general purpose integration method. 
However, there are a number of items that require further 
investigation and research. 

In particular, the area of stability in the pre- 
sence of stepchanges has not been dealt with adequately. 
Practice has shown that large errors can result if order 
q formulae are used and the stepsize is not held constant 
for at least q steps after a change. 

Because Nordsieck's method is based on polynomial 
approximations it appears that it is related to the 
problem of constructing stable polynomial filters. 

Hence, ideas from sequential smoothing and filter theory 


may lead to improvements. 


sabhibtioiisd “youse's Moletainss ce 
emeldoxd suisy’ isigint vileo ionic patylos 102 abodiom — ke , 

mt Ssgolevob ak borktom | ont | -anoissup3 rexsaoustRAs ak ol } 

nix bonimexs ai prtprstio astageta sit IT xssqeao i wis r 

to mofdsinomefqmt ns to Snilsvo teitd A «IIT tesqen> 7 


gest of¥ bas VI tetqsid ai nevip ef abodstom 2 ‘AseLebion © 
Jtoqque ativest seed? .V rstqedD al bstrogez sts attvaos 
eovisaeb yietinilsb borjem a'xosiebro# ssd+ motewlonoo edz 
-borssm aolssa1psini-seoaqivg isismep 6 as noitsrebrlenod F 
+oriduy? etbupex dadt amesi Fo 1ssdmun 5 sie ered? ,19svewoH 
.fdoxsetex bas noitspiteevak 
-syq ord oi ytiltdsde to Bex15 sds ,teluoidisg al 
-Yletsupsbs ditw tiseb need ton esd 2apnsiogede Yo eones 
rsbio ti tluasx ms a1o1rs epxzs6l jed3 awone asf eoisosi1t 
gastenoo blod gon. ei Ssteqete ods Das foay sis sailumxr02 p 
.Spnsd> 6 xstis eqste p sesel +8 208 wi 
fsimonyiogq no bsesd ef boijonm e'xoetabioll sausoed 


eit ot betslox 2k ti tends expeqqs ti enolteamixosqgs 


-exe3 LL? Inimonylog oldese paisouzsenop to mee . 
yxoonst pecia? bis paidsdooma isttnsypsa mor? esebs .90neH - 
.2insmevoigni oF Bsel beget Donte 


fi i 


_ a a! > ; i 


96 


Although Nordsieck's method handles first order 
systems of linear equations it is easy to modify the 
method to handle higher order equations directly. 
Superficially, there is some indication that solving 
high order equations directly may be faster than solv- 
ing the equivalent linear system. 

InNerve 1S also the question of how to handle 
Singularities. Clearly, if there is a singularity in 
the solution then polynomial methods are of little use, 
except for certain types of mild singularities such as 
a common point or an envelope for a family of solutions. 
Although, our experience has shown that Nordsieck's 
method integrates across such mild singularities it 
results in serious degredation of efficiency. Hence, 
examining special techniques to detect and deal with 
Singularities may be rewarding. However, if such tech- 
niques are included in a general purpose program then 
any advantage that may be gained is partially offset by 
additional overhead costs which are already substantial. 

As mentioned before our implementation of 
Nordsieck's method expends considerable effort in 
starting the solution. This problem is partly related 
LO, tie Stepchanging GLrors and lo, the accuracy of error 
estimates. Hence, any technique which decreases the 
step-changing errors produces more accurate error esti- 


mates and should improve the efficiency. 
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Methoden zur Integration gewoehnlicher Differen- 
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APPENDIX 


Program Listings 


The appendix contains the program listing of the 


Differential Equation Solver NDSK. It is followed by an 


example of a calling program for NDSK. 


Od 


$SIGN SEKE PRINT=PN FORM=BK 

ON AT 142312e45 ON TUE MAY 29/73 LAST ON AT 12:08 e49 
SSET LINECNT=50 

S$COPY *SOURCE* 


(axe) 


ANADDA 


NOAAAANDANA 


ANNDANADANANANAAANADANDANDIADDA 


DIFFERENTIAL EQUATION SOLVER NSOK 


THE FOLLOWING PROGRAM IS A COLLECTION GF SUBROUTINES 
WHICH SOLVE THE INITIAL VALUE PROBLEM IN THE FORM 


DY(T) 47 OT = F(TSY1eV¥2sceesYN) (1) 
GIVEN YIQCTO)J=YOs T=lelsecesNe 


SUBROUTINE NOSK(TOsTFeT POs HMAXs HMINs YOsNQGsECQs 
*ITODR+KDEQ,CQ) 


IN GENERAL» THIS SUBROUTINE ORGANIZES AND SUPERVISES THE 
INTEGRATION PROCESSe FIRST» THE INITIAL VALUES ARE TRANS— 
FERED INTO A COMMON BLOCKe NEXT» SOME PARAMETERS ARE 
CHECKED FOR FAULTSe THENs THE OTHER SUBROUTINES ARE 
CALLED IN THE REQUIREO ORDER TO CARRY OUT THE INTEGRATION 
PROCESSe FINALLYs THE TERMINAL VALUES OF THE SOLUTION 

ARE PLACED IN THE ARRAY YOo 


DOUBLE PRECISION A(10218)sYO0(10 )2TOs TFsTPOsHMAXeHMINe 
*TsTQsHeHWeTE 

DIMENSION Ca(3) 

COMMGN /BLKI/ZA 

COMMGN /BLK2/EC.ECTsI CLO eMAXRsNEQGQ ep KDELsEX!] eEX22EX3 
THE PARAMETERS HAVE THE FOLLOWING SIGNIFICANCESs 


TO IS THE INITIAL VALUE OF THE INDEPENDENT VARIABLE e 
TF IS THE FINAL VALUE OF THE INDEPENDENT VARIABLE’ 
HMAX IS THE MAXIMUM VALUE TO WHICH THE STEPSIZE MAY 

BE INCREASEDe ON DEFAULT» HMAX IS SET TO 106 
HMIN IS THE MINIMUM VALUE TO WHICH THE STEPSIZE MAY 

BE DECREASEDe ON DEFAULTs HMIN IS SET TO 0610-156 
TPO IS THE LENGTH GF THE INTERVAL BETWEEN 

PRINT—-POINTS. 
YO IS A VECTOR CONTAINING THE INITIAL VALUES OF THE 


DEPENDENT VARIABLEe WHEN CCNTROL RETURNS TO THE 
CALLING PROGRAM YO CONTAINS THE FINAL VALUE OF 
THE INDEPENDENT VARIABLEe 


NQQ IS THE NUMBER OF EQUATIONS IN THE SYSTEM (€ 1 )e 

ECQ IS THE ERROR TOLERANCE SUPPLIED BY THE CALLING 
PROGRAM. 

IO0R IS THE MAXIMUM ORDER TO BE USEDe THIS VALUE MUST 
BE GREATER THAN 1 AND LESS THAN OR EQUAL TO lé@e 

KDEQ IS THE MINIMUM NUMBER OF STEPS THAT THE PROGRAM 
TAKES BETWEEN ORDER CHANGES. 

CQ IS A VECTOR OF WEIGTHS USED FOR SCALING THE 
STEPSIZEe 
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) IS THE WORKING ARRAY. 
T=1s2e3e0ec0ee212sTS USED TO HOLD YC(T) AND ITS FIRST 
ELEVEN DERIVATIVES. 
) CONTAINS THE CORRECTION FOR A SINGLE ITERATIONe 
) CONTAINS THE ACCUMULATED VALUE OF THE CORRECTION. 
Y CONTAINS THE CORRECTION APPLIED AT THE 
PREVIOUS STEP. 
» CONTAINS THE BACKWARD OIFFERENCE OF THE 
CORRECT IONe 
» CONTAINS VARIQUS PARAMETERS 
) CONTAINS THE INITIAL VALUES OF YC T)e 


CHECK PARAMETER LIST FOR FAULTS» IF SUCH ARE FOUND» 
ATTEMPT TO CORRECT THESE» GR PRINT A MESSAGE AND RETURN 
CONTROL TO THE CALLING PROGRAMe ALSOse PLACE THE PARA- 
METERS IN THE WORKING ARRAY Ae 


1 ECT=ECQ 
EC=ECT 
MA XR=12 
NE @=NQQ 
KDEL=KDOEQ 
0G 40 I[=1.s3 
IF (CQCT) eL Te (003) ORC CQ(I) eGTe( 3e0)) CQCI)=140 
40 AC7+1.17)=OBLE(CQC(I)) 
DO 42 N=1-2NEQ 
42 A(N218)=YO(N) 
IF CHMINeLE 000) HMIN=1-0E-15 
ITF CHMAXe«LTeHMIN) HMAX=1.0 
CHMAX eGTeTPO) HMAX=TPO 


217)=TPO 
4517)=HMAX 
A(5217)=HMIN 
IF (CNEQ@eGEe le ANDOeNEQeLEC190) GO TO 50 
CALL RPRNT(0-O0D OO) 
RETURN 


IF 

AQ1 
A(2917)=TF 
AC3 

A 


SO IF( IOOReGE el eANDeITOOReLE*12) MAXR=IODR 


CALL SUBROUTINE START TO INITIALIZE THE WORKING ARRAY FOR 
EULER*S METHODe AND CALL SUBROUTINE PRNT TO PRINT THE 
INITIAL VALUES. 


T=TO 
CALL START (IOR+HsKSTEP»T) 
TQ=TO+DSIGN( TPO.H) 

60 CALL PRNT(0s0sIORsHs0.0D 007) 
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Cc INTEGRATE ONE STEP FORWARD WITH STEPSIZE Hs CALL 
Cc SUBROUTINE TEST TO DETERMINE THE PROPER STEPSIZE HWe 


100 CALL NORD (1sHsTsKSTEP.s. IOR) 
103 CALL TEST(KSTEPsT »He HW. IOR) 
A(6e17)=H 


IF HWs AS CALCULATEO BY TEST» IS GREATER THAN OR EQUAL 
TO H THE STEP IS ACCEPTED; OTHERWISE. THE STEP IS 
REPEATED WITH STEPSIZE HWe 


106 IFC CDABS(H)) eLEe(OABS(HW))) GO TO 120 
H=HW 


Cc 
Cc IF HW IS LESS THAN HMIN THE INTEGRATION IS TERMINATEDe 
Cc 


ANNADA 


IF( (DABS(HW))eGEeHMIN) GO TO 100 
108 TQ=T+tA(63:17) 

CALL SPRNT(CKSTEP.TO) 

GG TO is9 
120 KSTEP=KSTEP+41 

T=Tt+H 

H=HW 


IF THE NEXT STEP PASSES THE ENDPOINT BRANCH TO 
STATEMENT 1506 


TE=DABS(T-TF ) 
IF( TEeLEe(DABS(H))) GO TO 150 


IF THE NEXT STEP PASSES A PRINT-POGINT CALL PRNT TO 
INTERPOLATE ANDO TO PRINT THE SOLUTION AT THAT POINT. 


foverele) 


exerere) 


125 TE=DABS(T-TQ) 
IF( TE-DABS(H)) 13021353100 
130 CALL APRNT(KSTEPsHeTEsT) 
GO TO 140 
135 CALL NPRNT(K STEP. T) 
140 TQ=TQ+4DSIGN(TPO> H) 
GO TO 100 


THE FINAL VALUE OF Y(T) IS CALCULATED» PRINTEOs AND 
PLACED IN YOe 


150 IF(TFeLTeTO) TE=-TE 
CALL NORD(—-1 «TEsTeKSTEPsIOR) 
TO=Tt+TE 
CALL NPRNT(KSTEPsTO) 

159 CO 160 I=1.NEQ 

160 YOCI)=ACIe1) 
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999 RETURN 


END 


SUBROUTINE NORDCICALL»s HNs Ts KSTEPsNEWO) 


THE SUBROUTINE NORD INTEGRATES THE SYSTEM C47 
OF ORDINARY DIFFERENTIAL EQUATIONS EXACTLY ONE STEP 
FORWARD AT EACH CALL USING NORDSIECK*S METHOD WITH THE 
COEFFICIENTS DEFINED IN SUBROUTINE COEFTSe 

THE STEPSIZE HN MAY BE OIFFERENT FOR EACH CALLe THE 
QROER NEWO CAN VARY FROM 2 THROUGH 12% BUT SHOULD NOT 
CHANGE BY MORE THAN 1 FOR SUCCESSIVE CALLSe 

THE SUBROUTINE REQUIRES A DORIVING FUNCTION 


FEVAL (NEQoT eVe2F2TOeTFsHNsEC) 


TO BE SUPPLIED BY THE USER TO EVALUATE THE FUNCTION 
FOF ® (7 1m). 


DOUBLE PRECISION A( 10018) 0C€(3912) sASAVE (10216) eD(10) » 
*AE(102017)2Y(10)2F(10)sHLDs HSAVEs THe T sHN»sRHOeRHeTOeTF 
DIMENSION CP(12) 

EQUIVALENCE (A(6017) eHLD) os (AC 7217) sHSAVE) 9 (AC1217) TO) » 
¥CAC [S017 oT) oC ACL ol) oV (1) ) 

COMMON /BLKI/A 

COMMON /BLK2/ECeECT »s IOLDs MAXReNEGQe KDEL vEX1 sEX22EX3 


COMMON /BLK3/ASAVE2D 
COMMON /BLK6/C 


THE PARAMETERS HAVE THE FOLLOWING SIGNIFICANCE: 

ICALL IS A FLAGe IF ICALL IS A POSITIVE INTEGER 
Y(T4H) IS CALCULATED. IF ICALL EQUALS ZERO 
THE SOLUTION IS INTERPOLATED TO A PRINT— 
POINTe IF ICALL IS NEGATIVE THE SOLUTION IS 
EXTRAPOLATED TO THE THE END-POINT.- 


HN IS THE CURRENT STEPSIZE. 
a IS THE CURRENT VALUE OF THE INDEPENDENT VARIABLE’ 
NE WG IS THE OROER GF THE METHOC USED FOR THE 


PRESENT STEPe 


IOLD IS THE ORDER USED FOR THE PREVIOUS STEP. 

KSTEP IS THE STEP COUNT. 

KOOR IS THE VALUE OF KSTEP AT THE LAST ORDER CHANGE « 

ASAVE IS USED TO STORE A TO BE RETRIEVED IF THE STEP 
IS TO BE REPEATEDe 

AE IS USED TO STORE A WHEN THE SOLUTION IS TO BE 
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INTERPOLATED TO A PRINT-POINTe 


aN 


TH=T+HN 
NE=NEWO 
99 IFC ICALL) 103,100.102 


A IS SAVED TO BE REINSTATED WHEN THE INTERPOLATION TO A 
PRINT-POINT IS COMPLETED. 


160 DO 101 N=1,.NEQ 
DO 101 I=1.17 
101 AECNeTI=ACNeI) 
GO TO 131 
Cc CHECK WHETHER THE STEPSIZE OR THE ORDER IS TO BE CHANGED. 
Cc 


102 IFCNEWOeNE*TIOLD) CALL COEFTS(NEWO 202CeDeCP) 
IF (CABSCHLD)—-DABSCHN)) 10321034120 


Cc 
Cc A AND H ARE SAVED TO BE REINSTATED IF THE TESTS FAILe 
Cc 


ANNO 


103 HSAVE=HLD 

104 00 110 N=1s,NEQ 
ASAVE(Ne15 )=A(N2 14) 
ASAVE(N 016 )=A( N15) 
DO 110 I=1eNEWO 

110 ASAVE(Ns IJ=ACNS I) 
GO TO 126 


THE STEPSIZE IS TO BE REOUCED; HENCE» THE VALUES OBTAINED 
AT THE PREVIOUS STEP ARE REINSTATEDe 


ANNAN 


120 HLD=HSAVE 
DO 121 N=1-+NEQ 
A(N214)=ASAVE(N915) 
A(Ns 1S)=ASAVE (Ns 16) 
DO 121 I=1.NEWO 

121 A(NeI)=ASAVE(NeI) 


Cc 
Cc A CORRECTION IS CALCULATED IF H IS TO BE DECREASED. 
Cc 


IF(NEWOeLEe6) GO TO 126 
NW=NEWO+10 
CALL COEFTS(NWel eCeD0eCP) 
DO 124 N=1-2NEQ 
DO 124 J=3»eNEWO 
124 ACNeJ)=A(NoJ)tD( J—2) FASAVE (N913) 


Cc 
c THE BACKWARD OIFFERENCE OF THE HIGHEST ORDER SCALED 
Cc DERIVATIVE IS APPENDED TO CONTINUE THE INTEGRATION 
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Cc 


Cc 
126 


‘eke leke) ANAND 


ANNDAD 


128 
130 


106 


WITH HIGHER ORDER. 


IF CNEWOeLEeIOLD) GO TO 130 

DO 128 N=1.NEQ 
ACNeNEWO)=A( Ne 14)/C( 2eNEWO-1) 
IF(HNeEQGeHLD) GO TO 216 


THE VALUES IN THE WORKING ARRAY A ARE SCALED FOR A 


131 


132 


216 
218 


220 


270 


a 


276 


280 


282 
285 


STEPSIZE CHANGE. 


RH=0e10D O1 
RHO=HN/HLD 

DG 132 Y=2sNE 
RH=RHXx RHO 

O00 132 N=1.NEQ 
A( Noe I )=RH*XACNs I) 


THE PREDICTION IS DONE BY EFFECTIVELY MULTIPLYING THE 
VECTOR OF SCALED DERIVATIVES BY A PASCAL MATRIXe 


IF(NEeLEe1) GO TO 270 
DO 220 N=1;sNEQ 

DO 220 I=2>sNE 

DO 220 J=I .NE 
K=NE+4+I-J-—1 

ACN sKI=ACNoK) FAC Ne KH1 ) 


THE SOLUTION IS CORRECTEDe CORRECTIONS ARE CALCULATED 
UNTILL THE RESIDUAL IS LESS THAN THE ERROR TOLERANCE eo 
AT MOST THREE CGRRECTIONS ARE DONE PER STEPe 


DG 275 N=1>,NEQ 

AC Ne 1S)=ACNe 14) 

A(N214)=0-00 OO 

A(Ne13)=0.0D OO 

DO 285 L=13 

RH=HN 

CALL FEVAL (THs YeFsTOsTFsHNeoEC) 
IF (RHeNE-HN) GO TO 120 

OG 280 N=1.NEQ 


A(Ns13)=HN*XF (NI—AC(Ne 2) 
YON) =VON)4C0151)*ACN013) 
A(Ns2)=A(Ns2)+A(N913) 
AC(Ns14)=A(Ns 14) 4+A0N2 13) 


NQ=NEQ 

DO 282 N=1+NEQ 
IF(DABSCA(N213))eLT-DBLECEC/FLOAT(NEQ))) NQ=NQ-1 
CONTINUE 

IF(NQeLE°O0) GO TO 286 

CONTINUE 
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Cc 
Cc 
Cc 


AAA ANN|AADA 


ANANDA 


290 AC Ns16)=A(N0 14) 
294 ITOLD=NEWO 
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286 IFCICALL.EQe0) GO TO 296 


THE SCALED DERIVATIVES ARE CORRECTED. 


IFCNEWOeLTe3) GO TO 288 
DO 287 I=3»NEWO 
DO 287 N=1.NEQ 


287 ACNeTV=ACNoI D+C(1 oI) €ACNO1 4) 


THE BACKWARD DOIFFERENCE OF THE CORRECTION TERM IS 
CALCULATED TO BE USED TO ESTIMATE THE STEPSIZE FOR 
THE NEXT HIGHER ORDER. 


288 0G 290 N=1 »NEQ 


ASAVE(Ne 13 )=A(N2 14) 
AC(No14)=C(1sNEWO) *ACN 
—-A( Noel 


RE TURN 
THE SOLUTION AT THE PRINT-POINT IS PUT INTO A(N2i8}6 


296 DO 298 N=1.NEQ 


A(Ns18)=ACNs 1) 
DO 298 I=1.17 


298 AC(NsI)=AE( NOT) 
999 RETURN 


END 


SUBROUTINE TEST(KSTEPsT eHLDeHW» IGOR) 
SUBROUTINE TEST CHECKS THE STEPSIZE AND THE ORDERe 


OGUBLE PRECISION A( 10018) 2C(3012) eWTS(10) os TeHi DeHWe 
*RHO»sT1sT2sT3 

OIMENSION CP(12) 

COMMON /BLKIZ/A 

COMMON /BLK2/ECs ECT» IOLDeMAXReNEQG oe KDEL 0EX1 sEX22EXS 
COMMON /BLK4/WTS »CP.KODR eaKHC 

COMMON /BLK6/C 


THE FARAMETERS HAVE THE FOLLOWING SIGNIFICANCE: 
HW IS THE NEW STEPSIZEe 
HLD IS THE STEPSIZE WHICH HAS BEEN USED. 
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Cc 
Cc 
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ANN0 


60 
70 


80 


130 


137 


140 


200 


205 
206 


210 
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CHECK WHETHER THE ERROR CONSTANT HAS BEEN CHANGED BY 
FEVALs IF SOs THE NEW ERROR CONSTANT IS PLACED IN CPe 


THE MAXIMUM VALUE OF THE SOLUTION IS PLACED IN WTSe 


HW=HLD 

DG 70 N=1.NEQ 
WTS(N)=DMAX1 (WTSC(N) eDABSCACNol) )) 
IFC ECTeEQeEC) GO TO 130 

DO 80 I=1,12 

CP (T)=CP(I)*(CECT/EC) **2 

ECT=EC 


A CHECK IS DONE TO OETERMINE WHETHER THE STEPSIZE HAS 

BEEN REDUCED MORE THAN FOUR TIMES ON THE PRESENT STEP. 
KQ CONTAINS THE COUNT. IF SO.» THE ERRORS IN THE STGRED 
QUANTITIES ARE ASSUMED TO HAVE EXCEEDED THE ACCEPTABLE 
LEVEL ANDO THE INTEGRATION IS RESTARTED WITH EULER®*S 


METHODe 

IFC KSTEPeNE eKHCeOReKSTEPeLEe2) GO TO 140 
KQ=KQ+1 

IF(KQeLEe4eOReIOReLE2) GO TQ 200 

IOR=2 


CALL OPRNT(KSTEPs,IOReT) 

CALL RSTART CHW ) 

IF (DABS(HW)-eGEeDABSCHLDO)) HW=0e95*HLD 
RETURN 

KQ=0 


THE QUANTITIES REQUIRED TO CALCULATE HW ARE ACCUMULATED. 
T2 IS A MULTIPLIER SUCH THAT H/T2 IS THE REQUIRED HWe 


RHO=0-0 

DO 205 N=1 »NEQG 

RHO=RHOC+(A(N014) /WTSCN) )*¥2 

T2=A(9 017) *( CPC IGOR) *RHO4C( 3s IOR—1 ) #*¥ 2) ¥¥E XS) 


THE ORDER IS CHECKED ONLY EVERY KOEL STEPS. 

IF ( KSTEP—KODOR-—KDEL) 21022102220 

IF (CT2eGTe(0e8D 00) eANDeT2eL Te( 04100 01)) RETURN 
GO TO 300 


THE QUANTITIES T1 AND T3s REQUIRED TO FIND HW FOR THE 
NEXT LOWER OR HIGHER ORDERe ARE ACCUMULATED & 


220 T1=0e0D 00 


T3=0-0D O00 
KODR=KSTEP 
DO 240 N=1.NEQ 
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240 


ANNAA 


250 


Cc 

Cc 

S 
300 


305 
345 


NADA ANADDAD 


iS) 


T1=T1+C ACN» TOR) SWTSCN) )*¥*¥2 

T3=T3+(A(Ns 16)/WTS(N )) *%2 

RHO=0e1D 15 

IF (CIOReGTe2) RHO=A( 8517) *(CPCIOR—1)*T1*C(3esTOR—-1 )**2) 
*XXE XI 

T1=RHO 

RHO=02e1D 15 

IFC TOReLTeMAXR) RHO=A(10017 )*(CP( IOR+t1 ) #T3*C( 3. 1GR-1) 
¥EKO ) KKEXZ 

T3=RHO 

RHO=DMINI(T1.T22T3) 


RHO IS SUCH THAT H/RHO IS THE: REQUIRED STEPSIZE. 
IF RHO EQUALS T1 THE GRDER IS LOWERED. IF RHO EQUALS T3 
THE ORDER IS INCREASED. 


IF (RHOe GEe(0461D 01)) GO TO 300 

IF (DABS(T2—-RHO)—-00610*T2) 21092102250 
T2=RHO 

IF(T2eEQeT1) IGR=IOR-1 

IF(T2ceEQeT3) IOR=IOR1 
TOR=MAXC (IOR of ) 
EX1=14-0/(22¢0*FLGAT(IOR) ) 
EX2=1¢-¢0/(2e0*(1-e04FLOATCIOR))) 
EX3=120/(220*(2-e¢04FLOATCIOR))) 


THE NEW STEPSIZE HW IS CALCULATED. 


HW=0e9D OO*HLD/DMAX1(T22001D 00) 
HW=DSIGN(DMIN1(A(0 4017) 2 DABSCHW) )» HW) 
IF (DABS(HW) eLTeOABS(HLD)) KHC=KST EP 
RETURN 

END 


SUBROUTINE STARTC(CIORsHeKSTEPsT) 


SUBROUTINE START INITIALIZES THE WORKING ARRAY. 
THE MEANING OF THE ARGUMENTS IS THE SAME AS BEFORE: 


DOUBLE PRECISION A(10018) sC(3012)eASAVE( 10216) sD(10) >» 
*Y( 10) sFC10) oWTS(10) eHe TOs TF oT PO sHMAX eHMINaHLDs T 

DIMENSION CP (12) 

EQUIVALENCE (A(191)2Y(1)) 

EQUIVALENCE (A(1917) TO) o(AC2017) s TF) 0 (AC4 017) sHMAX) 2 
*(A(6e017)2HLD) 


; 7 ! : : 


: 7 es ie a: 
Pe ioe Ba 
‘ih ei sei. 
(Sek(1-ADLeEI DEL T#CE-AOTIGIIS bi 


( P=AO] -EDI¥E TAC TVROL GD INIT CORIAZOHF (FR , 


Ey esau One ai "a38 2 on 


| a 
(tH 


nO Laeag aed aoe 
en Sy. 


((tHOF 


9 GkeOeST)IXAMGY GIF 
CH of CMH) 2BAG (Teg) oni 
QSTSASOHM (I CIHI2SAGeT 


(T VATS etieROLITAAT2 BME TwaRsaVe 
YAASA DMIAAOW SHT 2SRIAATTIMI TAAT2 SMES 

«390836 2A AMA2 BHT 21 Senn: GORA SHT 90 wathate nT 
(OL30 ce 4OLeOLPBVAZAs( ST eE3De (8S OLJA VWOtaI >a 
: ieee UROMNEEN RIE: f 


nieuekcre 
UKAMHE (Te DADS EAT eC te SPAR MOPEET ES EDAD, 


ANNAN 


Cc 
e 
Cc 


A ANDANDN 


110 


COMMON /BLKIZA 

COMMON /BLK2/ECs, ECT» IOLDeMAXRs NEQe KDEL »s EX1 es EX2sEX3 
COMMON /BLK3/ASAVE 2D 

COMMON /BLK4/WTS CPs» KOORsKHC 

COMMON /BLK6E/C 


IOR=2 

KSTEP=0 

DO 102 I=1,.10 
O(1)=0-00 00 


102 WTS(I)=0e1D 01 


THE INITIAL VALUESe AS PROVIDED BY THE CALLING PROGRAM, 
ARE PLACED IN THE WORKING ARRAY. 


DO 105 N=1+NEQ 
ACNe1)=A(N018) 


16S ASAVE(Ns1)=AC Net) 


110 


THE ERROR CONSTANT TO BE USEO INITIALLY IS CALCULATED. 


CALL CCEFTS(IGRs—-12CeDeCP) 
C(C2e1)=0-1D0 O1 
C(3s1)=C0.10 01 


C(1+*eI)=0-0D 00 
C(2eT)=001D O14C(2,I-1) 
CPCI )=CCPCI) SECT ) **2 
C(3sT)=C(3e1-1)*C(2sT) 


C(122)=021D0 O1 


THE REQUIRED PREPARATIONS ARE MADE HERE TO START THE 


INTEGRATION PROCESS WITH THE CLASSICAL EULER METHOD 


USING THE PRESENT VALUES OF THE VARIABLESs DEPENDENT 
AND INDEPENDENTs AS THE INITIAL VALUESe 


ENTRY RSTART(H) 


1QLD=2 

KODR=KSTEP 
KHC=KSTEP 
H=A(5917)**00e25 
IF ( TFelTetTO) H=—-H 
HL D=H 

A(7217)=H 

CALL FEVAL(T sYoF 2s TOoTF sHsEC) 
DO 275 N=1.NEQ 
ACN e2) =H*F (N) 

DO 275 I=3212 
A(NsI)=C-0D O00 
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275 ASAVE(Ns1I)=0-0D OO 
ExX1=0e25 
EX2=1-.0/6.0 
ExX3=0e125 
CALL COEFTS(IOR.02C.DsCP) 
RETURN 
END 


SUBROUTINE COEFTS(CIOORs IXeCsDsCP) 

SUBRGUTINE COEFTS ASSIGNS THE PROPER COEFFICIENTS TO BE 
USED FOR THE IODR*TH ORDER METHOD. 

DOUBLE PRECISION C(3012)3D(10) 

DIMENSION CP(12) 


IFC IX) 300310024200 


100 GO TO (102:359425260728292100119312) »sI0DR 

1 C(121)=0-1D O1 
RETURN 

2 C(121)J=0e5D 00 
RETURN 

3 C(1201)=0e5D 00 
C(123)=0e5D 00 
RETURN 

4 C(1.1)=0-4166666666666667D 00 
C(123)=0e75D O00 
C(1 44) =0-1666666666666667D 00 
RE TURN 

§ C(121)=0e375D 00 
C( 1.463) =0-9166666666666667D 00 
C(144)=0-33333333333333330 00a 
C(155)=0 -4166666666666667D-—-01 
RETURN 

6 C(1.61)=0-3486111111111111D 00 
C(1+e3)=0-1041 666666666667D Ol 
C(1.4)=00e4861111111111111D OO 
C(1 25)=0.1041666666666667D 00 
C( 1.66) =0-e8333333333333333D-02 
RETURN 

7 C€1.1)=0232986111111111110 CO 
C(143)=0-1141 6666666666670 Ol 
C(124)=0e625D 00 
C(1.5)=00e17708333333333330 co 
C(1.6)=0-250-01 
C(107)=0.1388888888888889D-0 
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3155919312169312D 00 
1225D 01 

7518518518518518D 00 
25520833333333330 00 
©48611111111111110-01 
e48611111111111110-02 
19841269841 26984D- 03 


030422453703 70370D 40 
¢©12964285714285710 Ci 
©86851 851851851850 Co 
335763888 8888889D 00 
e777 7T77777777780-01 
e1064814814814815D-01 
e 793650793650 7936D-03 
¢e24801587301587200— 04 


©29486800044091 71D 00 
e13589285714285710 01 
© 9765542328042328D 00 
e4171875D 00 

01113541 6666666670 00 
e1875D0-01 

© 1934523809523809D-02 
e1116071428571428D-03 
062755 7319223985 80-05 
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=0 e2 869754464285714D 00 

=02e1414484126984126D Ql 
c(1 =0-1077215608465608D Ol 
C(125)=02498567 0194003505D 00 
C(C1s6)=00e1484375D 00 
C(197)=0e29060570987654270-01 
C(1+8)=0-e37202380952380 950-02 
C(109)=0 -29968 584656 08465D— 03 
C(1 910)=0-1377865961199294D-04 
C€1.11)=062 7557319223985 890-06 
RETURN 


C(€121)=0022801895964439367D 00 
C(193)=00e1464484126984123D O1 
CO154)=0211715145502645460 O01 
C(125)=0-5793581900352713D CO 
C(1+6)=0-1883228615520275D 00 
C107 )=0 041430 362654320920-01 
C(128)=0.62111441 798941 740-02 
C0129 )=0 e62520 66798941 796D-—-03 
C(1910)=0-4041 740152851263D0-04 
C(1+11)=0. 151565255 73192240-05 
C(1012)=0-25052108385441 750-07 
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RETURN 


GO TO (16017218e19220%s21222)sI0DR 


D(1)=0-68247126436781650-01 
O0( 2)=04e7782567049808428D-01 
0(3)=0-4992816091954023D-01 
0(4 )=0-1709770114942529D-01 
RETURN 

0(€1 )=0-4841597796143251D-01 
0(2)=0 .67148760330578510-01 
0(3)=0-56387741046831950-01 
D(4)=0 e 2871900826446 2810-01 
D(S5)=0 282529 843893438025D-02 
RETURN 

D0(1)=0.32694892473118280-01 
0(2)=0-53419952210274790-01 
0(3)=0-5542254704301075D0D-01 
0(4 )=0 3743839605734 7670-01 
5)=0216109431003584230-01 
Y=0 -4031458013311852D-02 
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e41008690828402370-01 
e 4131030161078238D-01 
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e 1479099159387621D0-01 
0 « 4930856403212172D0-02 
0-9861905762012599D0-03 
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)=0-202 48541 8566051 70-01 
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THE COEFFICIENTS OF THE LOCAL TRUNCATION ERROR TERMS ARE 
PLACED IN CPe CP(P) WILL CONTAIN THE COEFFICIENT FOR THE 
P*TH ORDER METHOD. 


300 CP(1)=0-5 

CP( 2)=0- 08334 
CP (3)=0 -04167 
CP(4)=0.02639 
cP(S5)=0-01875 
CP(6)=0-01427 
CP(7)=0-01137 
CP (8)=054« 00936 
CP(9}=0.00789 
cP(10)=0.00678 
CP (11)=0 00527 
cCP(12)=0.00412 
RETURN 


SUBROUTINE PRNTCIPXsKSTEPs IOQRDEReHeHP oT) 


SUBROUTINE PRNT IS A COLLECTION OF ALL WRITE AND THEIR 
ASSOCIATED FORMAT STATEMENTSe TO AVOID MULTIPLE TESTS» 
SEPARATE ENTRY POINTS ARE PROVIDED FOR MOST WRITE STATE- 
MENTSe THE EXPLANATIONS ARE OMITTEOs THESE CAN BE 
CBTAINED BY CHECKING THE ASSOCIATED FORMAT STATEMENTS. 
IPX IS A FLAGe IF IPX=0 YeTs AND KSTEP ARE PRINTEDe 

IF IPX=1 THE SOLUTION IS INTERPOLATED TO A 

PRINT—POINTe IF IPX=-1 THE INTEGRATION TERMINATESe 
HP IS THE OIFFERENCE BETWEEN THE CURRENT VALUE OF 

T AND THE NEXT PRINT-POINTe 


DOUBLE PRECISION AC 10018) +Y(10) sF (10) sYP(10) oT sHeHPs TZ 
COMMON /BLKIZA 

COMMON /BLK2/ECsECTs IOLDsMAXRs NEQ eKDEL » EX1 sEX22EX3 
EQUIVALENCE (A(191)sY)2CAC1 02) 0F) 9(AC 1218) oYP) 
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950 FORMAT (* %¢1693XsF1200¢503XeS5(D19)e L201 XIS27X2S(D196e1291X) ) 

951 FORMAT(*® INTEGRATIGN TERMINATEDeH HAS BECOME TO SMALL®*) 

952 FORMAT(* PARAMTERS ARE IMPROPERLY SPECIFIED.INTE— 
*GRATION NOT STARTED ‘*) 

953 FORMATC* OTHE GUTPUT COLUMNS ARE *//* «STEP T YC1) 
*Y¥(2) eco YCNEQ) */77°O0 *¢) 

954 FORMAT (°* AT START T = %09F12¢5910Xe*STEPSIZE = %s 
*¥F12¢525Xe*GROER = *.I16//) 

955 FORMAT #0 es 30Xs*AT T = *%sE12e5e5Xe *THE ORDER CHANGED 
*TO %eTGe2SXe*KSTEP = *,16//) 

956 FORMAT( 70% s30Xs"AT T = %sE120e595Xe*H CHANGED FROM ‘%-¢ 
*E1L120¢725Xs*TO ®eEL2ce7sS5Xe*KSTEP = *%,16/7) 

957 FORMAT(*O THE INTEGRATION IS RESTARTED USING EULER'S 
*METHOD # OF STEPS %e16e* ORDER *3160° T *%sD0200e13s° 
* H *,D200e13// 


NQ=NEQ 

IFCKSTEPeGTe0) GO TO 4 

WRI TE( 6,953) 

WRITE(62954) TseHs IGRDER 
4 IFCIPX) 20230210 


ENTRY RPRNT(T) 
S WRITE(6 2952) 
RETURN 


ENTRY APRNT(KSTEPsHs HP eT) 

10 IFC HeLTe(00e0D 00)) HP=—DABS(HP) 
IORDER=I0LD 
CALL NORD(0sHP»s TsKSTEPs IORDER) 
TZ=HP+T 
WRITE(62950) KSTEPsTZs(YP(I)sI=1sNQ) 
RETURN 


ENTRY OPRNTC(KSTEP>s IORDER oT ) 
12 WRITE(62955) TsIORDERsKSTEP 
RETURN 


ENTRY HPRNT(KSTEPeHsHP oT) 
15 WRITE(62956) TeHeHP»KSTEP 
RETURN 


ENTRY TPRNT(KSTEPsIORDEReTsH) 
16 WRITE(6e957) KSTEPsTIORDEReT sH 
GO TO 30 


ENTRY SPRNTC(KSTEP oT) 
20 WRITE(62951) 


ENTRY NPRNTC(KSTEP oT ) 
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TZ=T 
WRITE(6956) 
RE TURN 

END 


KSTEPsTZs(YC(I) »eI=15NQ) 
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$S$COPY *SOURCE* 
SAMPLE PROGRAM TO CALL NDSK 


CALLING PROGRAM TO INTEGRATE: THE 
SINE —- COSINE OIFFERENTIAL EQUATION 


ANDANAD 


DOUBLE PRECISION YO(10)2TO.s TFs TPO sHMAXsHMIN 
DIMENSION CQ(3) 
COMMON /BLK/NOF «ECT 
905 FORMAT(*O THE # OF FUNCTION EVALUATIONS = %es1625XsE10-5) 
TO=0-e0D 00 
TF=04e2D 03 
TPO=0e1D 02 
HMAX=0e5D O01 
HMIN=0-10-12 
YO(C1)=0-20D 00 
YO(C2)=0.21D 01 
NEQ=2 
ECT=0e1E-06 
IODR=10 
KDEL=0 
CQ(1)=1.1 
CQ(2)=1e2 
CQ(3)=1-3 
NOF =O 
CALL NOSK( TOs TFs TPOs HMAX sHMINe YOs NEGs ECT oI GCORs KDEL » CQ) 
WRITE(62905) NOFs ECT 
STOP 
END 


SUBROUTINE FEVAL(CTsYsF+TO 
DOUBLE PRECISION F(10) YC 
COMMON /BLK/NOF>s ECT 
NOF=NOF +1 

FC1)=V(2) 

F(2)=-Y (1) 

RETURN 

END 


FeHNeEC) 
YeoTOo@lF sHNoTl 


(09 ¢ J90M ¢HOO le TID yOSMeOY eMIMHa XA 


T MH ESTIOTS COTY SGOT ii i 7 >: 


